Función de Derivadas parciales
partial_dev <- function(x,i,fun,h=0.01){
e <- x*0 # crea un vector de ceros de la misma longitud de x
e[i] <- h
y <- (fun(x+e)-fun(x-e))/(2*h)
return(y)
}
Para evaluar el gradiente de una función f en x hay que obtener cada una de las derivadas parciales de f en x. Con el código anterior y la función mapply() se puede lograr esto fácilmente, como se muestra a continuación:
num_grad <- function(x,fun,h=0.01){
# x: punto del espacio donde se debe evaluar el gradiente
# fun: función para la que se desea calcular el gradiente en x
# h: es el tamaño de ventana para el cálculo de la derivada numérica
d <- length(x)
y <- mapply(FUN=partial_dev,i=1:d,MoreArgs=list(x=x,h=h,fun=fun))
return(y)
}
Función de Calculo de \(\frac{\partial}{\partial x_{t}} \nabla f^{T}\)
deriv_grad <- function(x,fun,i=1,h=0.01){
# x: punto en el que se evalúa el gradiente
# fun: función para la cual se calcula la derivada del gradiente respecto a la íesima componente
# i: i-ésima componente del vector x con respecto a la que se deriva
e <- x*0 # crea un vector de ceros de la misma longitud de x
e[i] <- h
y <- (num_grad(x+e,fun=fun,h=h)-num_grad(x-e,fun=fun,h=h))/(2*h)
return(y)
}
Función para Calcular la matriz Hessiana
matriz_hessiana <- function(x,fun,h=0.01){
# x: punto en el que se evalúa la matriz hessiana
# fun: función a la que se le calcula la matriz hessiana en x
# h: es el tamaño de ventana para el cálculo de la derivada numérica
d <- length(x)
y <- mapply(FUN=deriv_grad,i=1:d,MoreArgs=list(x=x,h=h,fun=fun),SIMPLIFY = TRUE)
return(y)
}
optimizador_mult_numdev <- function(x0,fun,max_eval=100,h=0.01,eta=0.01){
x <- matrix(NA,ncol =length(x0), nrow = max_eval)
x[1,] <- x0
for (i in 2:max_eval){
num_grad_fun <- num_grad(x[i-1,],fun,h)
H <- matriz_hessiana(x[i-1,],fun,h)
cambio <- - eta*solve(H)%*%num_grad_fun
x[i,] <- x[i-1,] + cambio
cambio_opt <- sqrt(sum((x[i-1,]-x[i,])^2))
if (cambio_opt<0.00001){
break
}
}
return(x[1:i,])
}
f_rosenbrock2d <- function(x1,x2){
# Esta versión de la función se usa para graficar.
y <- 100*((x2-x1^2)^2)+(1-x1)^2
return(y)
}
f_rosenbrock2d_vec <- function(x){
# Versión vectorizada. Es la que se utiliza para optimizar.
x1 <- x[1]
x2 <- x[2]
y <- 100*((x2-x1^2)^2)+(1-x1)^2
return(y)
}
Proceso de Optimización con condición Inicial Aleatoria \((x_{1},x_{2})^T\) con \(x_{i} \sim U_{i} (-2.048,2.048),\hspace{1mm}i=1,2\)
Solución
set.seed(123456)
x1_ros_inic <- runif(1, min =-2.048, max = 2.048)
set.seed(789)
x2_ros_inic <- runif(1, min =-2.048, max = 2.048)
x0_rosenbrock <- c(x1_ros_inic,x2_ros_inic)
sol_ros <- optimizador_mult_numdev(f_rosenbrock2d_vec,x0=x0_rosenbrock,eta=1.3)
Graficamente la evolución de la solución:
Función para Generar las imagenes en cada iteración
library(stringr)
library(magick)
plot_generate_function <- function(out_dir,sol_func,name_function,x1,x2,z){
n_length <- 150
X <- expand.grid(x1,x2)
Z <- matrix(z,ncol=n_length,nrow = n_length)
tamaño_iter <- dim(sol_func)
n_iter <- tamaño_iter[1]
for(i in 2:n_iter){
png (paste0(out_dir,"/",str_pad(i, 3, pad = "0"), ".png"))
contour(x1, x2, Z, las=1,
xlab = expression(x[1]), ylab = expression(x[2]),
main = paste("Función de ",name_function ," iter:",i),
sub = "Curvas de nivel de la función")
lines(sol_func[1:i,], type="b",cex=1.5,col="red")
dev.off()
#file.show(paste0(out_dir,"/",i,".png"))
}
png (paste0(out_dir,"/",str_pad(n_iter+1, 3, pad = "0"), ".png"))
contour(x1, x2, Z, las=1,
xlab = expression(x[1]), ylab = expression(x[2]),
main = paste0("Función de ",name_function),
sub = "Curvas de nivel de la función")
lines(sol_func[1:n_iter,], type="b",cex=1.5,col="red")
points(sol_func[n_iter,1],sol_func[n_iter,2],pch=8,cex=2, type = "p")
}
Generar las Imagenes
plot_generate_function("~/Tarea 07/img_ros_grad",sol_func=sol_ros,"Rosenbrock",x1=x1_ros,x2=x2_ros,z=z_ros)
Función para crear el Gif Animado
save_to_gif <- function(dir,out_dir,out_name){
imgs <- list.files(dir, full.names = TRUE)
img_list <- lapply(imgs, image_read)
## join the images together
img_joined <- image_join(img_list)
## animate at 2 frames per second
img_animated <- image_animate(img_joined, fps = 2)
## view animated image
img_animated
## save to disk
image_write(image = img_animated,
path = paste0(out_dir,"/",out_name,".gif"))
}
Guardar el Gif Animado
save_to_gif("~/Tarea 07/img_ros_grad","~/Tarea 07","Rosenbrock_Gradiente")
Mostrar Gif
Se observa entonces que la otpitmización por descenso por gradiente con condición inicial aleatoria, alcanza su valor mínimo en las coordenas (1,1), es decir el algoritmo alcanza el mínimo global de la función.
Graficamente la Solución
Observemos el comportamiento de la optimización más de cerca:
Generar las Imagenes
plot_generate_function("~/Tarea 07/img_ras_grad",sol_func=sol_ras,"Rastrigin",x1=x1,x2=x2,z=z)
Generar Gif Animado
save_to_gif("~/Tarea 07/img_ras_grad","~/Tarea 07","Ranstrigin_Gradiente")
Mostrar Gif
A pesar de que estuvo muy cerca de encontrar el mínimo global en las coordenadas (0,0), el método quedo atrapado en un mínimo local, encontrando como optimo las coordenadas (0,2).
f_rosenbrock2d_vec_inv <- function(x){
# Versión vectorizada. Es la que se utiliza para optimizar.
x1 <- x[1]
x2 <- x[2]
y <- 100*((x2-x1^2)^2)+(1-x1)^2
return(-y)
}
Graficamente la Solución
Crear las imagenes para el Gif Animado
plot_generate_function_ga("~/Tarea 07/img_ros_ga",sol_func=opt_ros_ga,"Rosenbrock",x1_ros, x2_ros,Z_ros)
Generar Gif Animado
save_to_gif("~/Tarea 07/img_ros_ga","~/Tarea 07","Rosenbrock_ga")
Mostrar Gif
El algoritmo evolutivo encuentra perfectamente la solución en las coordenadas (1,1).
f_rastrigin2d_vec_inv <- function(x){
# Versión vectorizada. Es la que se utiliza para optimizar.
x1 <- x[1]
x2 <- x[2]
y <- 20+x1^2-10*cos(2*pi*x1)+x2^2-10*cos(2*pi*x2)
return(-y)
}
Gráficamente la Solución
contour(x1_ras, x2_ras, Z_ras, las=1,
xlab = expression(x[1]), ylab = expression(x[2]),
main = "Función de Rastrigin",
sub = "Curvas de nivel de la función")
lines(opt_ras_ga@population, type="b",cex=1.5,col="red")
Crear las imagenes para el Gif Animado
plot_generate_function_ga("~/Tarea 07/img_ras_ga",sol_func=opt_ras_ga,"Rastrigin",x1_ras, x2_ras,Z_ras)
Generar Gif Animado
save_to_gif("~/Tarea 07/img_ras_ga","~/Tarea 07","Ranstrigin_ga")
Mostrar Gif
Se observa que el algoritmo genético es más robusto a mínimo locales y alcanza su mínimo global con satisfacción.
#Función de optimización PSO
opt_pas_pso_ras <- psoptim(
par = c(NA,NA),
fn = f_rastrigin2d_vec,
lower = c(-5.12,-5.12),
upper=c(5.12,5.12),
control = list(trace.stats=TRUE, trace=1, maxit = 150)
)
## S=12, K=3, p=0.2297, w0=0.7213, w1=0.7213, c.p=1.193, c.g=1.193
## v.max=NA, d=14.48, vectorize=FALSE, hybrid=off
## It 10: fitness=2.186
## It 20: fitness=1.351
## It 30: fitness=0.5029
## It 40: fitness=0.2419
## It 50: fitness=0.07325
## It 60: fitness=0.07254
## It 70: fitness=0.008249
## It 80: fitness=0.0001511
## It 90: fitness=1.973e-05
## It 100: fitness=1.973e-05
## It 110: fitness=2.112e-06
## It 120: fitness=4.153e-07
## It 130: fitness=4.153e-07
## It 140: fitness=6.349e-08
## It 150: fitness=1.037e-08
## Maximal number of iterations reached
Valor Optimo:
opt_pas_pso_ras$par
## [1] 1.951811e-06 6.960436e-06
Generar las imagenes
plot_generate_pso("~/Tarea 07/img_ras_pso",sol_func=opt_pas_pso_ras$par,"Rastrigin",x1_ras, x2_ras,Z_ras,opt_pas_pso_ras)
Generar Gif Animado
save_to_gif("~/Tarea 07/img_ras_pso","~/Tarea 07","Ranstrigin_pso")
Mostrar Gif
El algoritmo PSO muestra como las posibles soluciones candidatas convergen al mínimo global, siendo este también un buen método de optimización bastante robusto a mínimo locales.
#Función de optimización PSO
opt_pas_pso_ros <- psoptim(
par = c(NA,NA),
fn = f_rosenbrock2d_vec,
lower = c(-5.12,-5.12),
upper=c(5.12,5.12),
control = list(trace.stats=TRUE, trace=1, maxit = 150)
)
## S=12, K=3, p=0.2297, w0=0.7213, w1=0.7213, c.p=1.193, c.g=1.193
## v.max=NA, d=14.48, vectorize=FALSE, hybrid=off
## It 10: fitness=0.3783
## It 20: fitness=0.02858
## It 30: fitness=0.003881
## It 40: fitness=2.138e-05
## It 50: fitness=2.138e-05
## It 60: fitness=2.138e-05
## It 70: fitness=4.776e-06
## It 80: fitness=2.491e-06
## It 90: fitness=2.491e-06
## It 100: fitness=1.403e-06
## It 110: fitness=1.281e-06
## It 120: fitness=1.278e-06
## It 130: fitness=1.278e-06
## It 140: fitness=1.278e-06
## It 150: fitness=1.278e-06
## Maximal number of iterations reached
Valor Optimo:
opt_pas_pso_ros$par
## [1] 1.001130 1.002262
Generar Gif Animado
save_to_gif("~/Tarea 07/img_ros_pso","~/Tarea 07","Rosenbrock_pso")
Mostrar Gif
Igualmente el algoritmo PSO muestra como las posibles soluciones candidatas convergen al mínimo global, siendo este también un buen método de optimización bastante robusto a mínimo locales.
Gráficamente la solución
contour(x1_ros,x2_ros,Z_ros)
lines(opt_ros01$member$pop,type="p",pch=2,col="red",lwd=3)
Generar las Imagenes para el Gif
plot_generate_function_ED("~/Tarea 07/img_ros_ed",sol_func=opt_ros01 ,"Rosenbrock",x1_ros, x2_ros,Z_ros)
Generar Gif Animado
save_to_gif("~/Tarea 07/img_ros_ed","~/Tarea 07","Rosenbrock_ed")
Gif Animado:
El método de evolución diferencial converge igualmente bien, encontrando el mínimo global con satisfacción.
Gráficamente la solución
contour(x1_ras, x2_ras, Z_ras)
lines(opt_ras01$member$pop[], type = "p", pch = 2, col ="red", lwd = 3)
plot_generate_function_ED("~/Tarea 07/img_ras_ed",sol_func=opt_ras01 ,"Rastrigin",x1_ras, x2_ras,Z_ras)
Generar Gif Animado
save_to_gif("~/Tarea 07/img_ras_ed","~/Tarea 07","Rastrigin_ed")
Gif Animado:
El Algoritmo cae en un mínimo local muy cerca del global en las coordenadas (1,0)